######################################################################
# Time-stamp: <2005-07-15 09:03:21 patrick>
#
# Package: EP::Core::Dataset::Data::seq
#
# Sequence data subclass for the Data object (to sit inside <EP::Common::DBI::Dataset>)
######################################################################
package EP::Core::Dataset::Data::seq;

use strict;
use base 'EP::Core::Dataset::Data';

use EP::Config;
use EP::Common::General;

use Bio::SeqIO;
use PDL;

# sub: initialize
#   Lazy-loading routine to instantiate the sequence data object associated with a certain dataset. If an existing 
#   sequence IO stream is passed as the argument, this will be used to create the sequence data object.
# parameters:
#   Can /optionally/ take one parameter, a hash table with the string "sequences" as the key and 
#   an object reference as the value, in which case the sequence data object will be created.
#   The passed reference should be to an instance of a subclass of <Bio::DB::GFF>.
# returns:
#   An instance of object type EP::Core::Dataset::Data::seq
sub initialize {
  my $self = shift;
  my %params = @_;

  if ( ref( $params{sequences} ) =~ /^BioMart::ResultTable/ or
       ref( $params{sequences} ) =~ /^Bio::DB::GFF/ ) {
    $self->{sequences} = $params{sequences};
  }

  return $self;
#  $self->process_metadata() if not $self->{nr_sequences};
}

# sub: process_metadata
#   Process sequence specific metadata (none yet at the moment).
sub process_metadata {
  my $self = shift;

  my $seqs = $self->sequences;

  my @L;
  my $s;
  push @L, $s->length() while $s = $seqs->next_seq();

  $self->{nr_sequences} = scalar @L;
  my $l = pdl @L;

  my $M = { nr_sequences => scalar @L,
	    min_length => min($l),
	    max_length => max($l),
	    avg_length => int(sclr(average($l)))
	  };

  foreach my $metadata_type_name ( keys %$M ) {
    my $metadata_type = EP::Common::DBI::Dataset_Metadata_Type->find_or_create ( { name => $metadata_type_name } );
    my $metadata = $self->{owner}->add_to_metadata ( { type => $metadata_type->id } );
    my $value = $M->{$metadata_type_name};
    $metadata->add_to_values ( { value => $value } );
  }

  return $M;
}

# sub: sequences
#   Returns a <Bio::SeqIO> object for access to sequences
sub sequences {
  my $self = shift;

  my $dataset_filename = $EP::Config::h->{EPNG_USRDATA_PATH} . "/" . $self->{owner}->folder->location . "/" . $self->{owner}->filename . ".fasta";
  return Bio::SeqIO->new ( -file   => $dataset_filename, -format => 'fasta' );
}

# sub: start_bp
#   Set/Retrieve the starting basepair position for the dataset object of type seq.
#   Needed when retrieving for instance the upstream region of certain sequences, instead of
#   the whole sequence.
# parameters:
#   Can /optionally/ take one parameter, the starting basepair position.
# returns:
#   The starting basepair for the dataset object of type seq
# see also:
#   <end_bp>
sub start_bp {
    my $self = shift;
    my $bp = shift;

    if ( $bp ) {
	die "Did not get a numerical value for start basepairs" if EP::Common::General::data_type( $bp ) != NUMERIC;
	$self->{start_bp} = $bp;
    }

    return $self->{start_bp} if exists $self->{start_bp};
}

# sub: end_bp
#   Set/Retrieve the ending basepair position for the dataset object of type seq.
# parameters:
#   Can /optionally/ take one parameter, the starting basepair position.
#   Needed when retrieving for instance the upstream region of certain sequences, instead of
#   the whole sequence.
# returns:
#   The ending basepair for the dataset object of type seq
# see also:
#   <start_bp>
sub end_bp {
    my $self = shift;
    my $bp = shift;
    
    if ( $bp ) {
	die "Did not get a numerical value for end basepairs" if EP::Common::General::data_type( $bp ) != NUMERIC;
	$self->{end_bp} = $bp;
    }

    return $self->{end_bp} if exists $self->{end_bp};
}

sub bp_length {
  my $self = shift;

  my $start_bp = $self->start_bp();
  my $end_bp = $self->end_bp();

  if ( not $start_bp ) {
    my $seqs = $self->sequences;
    my $l = 0;

    while ( my $seq = $seqs->next_seq ) {
      my $s = $seq->length;
      $l = $s if $l < $s;
    }

    return $l;
  }

  my $bp_length = $end_bp - $start_bp;

  if ( $bp_length < 0 ) {
    return $bp_length + 2;
  } else {
    return $bp_length + 1;
  }
}

1;
